Jamming at zero temperature, zero friction, and finite applied shear stress 
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Via molecular dynamics simulations, we unveil the hysteretic nature of the jamming transition of 
soft repulsive frictionless spheres, as it occurs varying the volume fraction or the shear stress. In a 
given range of control parameters the system may be found both in a flowing and in an jammed state, 
depending on the preparation protocol. The hysteresis is due to an underlying energy landscape 
with many minima, as explained by a simple model, and disappears in the presence of strong viscous 
forces and in the small a limit. In this limit, structural quantities are continuous at the transition, 
while the asymptotic values of two time quantities such as the self-intermediate scattering function 
are discontinuous, giving to the jamming transition a mixed first-order second-order character close 
to that found at the glass transition of thermal systems. 

PACS numbers: 45.70-n,83.50.Ax,83.10.Tv 
Keywords: 



Athermal systems such as foams and granular materi- 
als undergo a transition from a disordered fluid like state 
to an amorphous solid, known as jamming transition, as 
the density is increased and/or the applied shear stress 
decreased. The properties of this transition, which is of 
interest due to its connections to the glass transition of 
thermal systems, and for its relation to important geo- 
physical phenomena such as avalanches and earthquakes, 
have not been fully clarified. Focussing on a system of 
repulsive soft spheres, Ref. [l| showed that the jamming 
transition at zero temperature and zero applied shear 
stress a occurs at a well defined point of the volume frac- 
tion <j) H, identified with the random close packing vol- 
ume fraction rcp (point J). The jamming transition at 
point J has a mixed first-order second-order character, 
whose origin is partially explained by simple geometrical 
percolation-based models [3J. Some quantities, such as 
the pressure P and the shear modulus G, are zero be- 
low the transition, and continuously grow as power laws 
above the transition, while the mean contact number Z , 
which is also zero below the transition, discontinuously 
jumps to isostatic value Z = Zi SO at the transition, and 
then grows as a power law. Zi SO = 2D, where D is the 
dimensionality, is the theoretical minimum value of the 
mean number of contacts required for mechanical stabil- 
ity [j], 0, HI- Finally, the shear viscosity in the limit of 
small applied shear stress exhibits a divergence as the 
jamming point is approached [6j. It has not yet been 
fully explored how this scenario changes in the presence 
of a finite value of the shear stress [H, H, 0] • 

In this Letter, we investigate via molecular dynam- 
ics simulations the jamming transition at a finite value 
of the applied stress, i.e. in the <p-a plane, discovering 
that this is actually hysteretic. In a region of the 4>-a 
plane the system may be found both in a flowing and in 
a jammed state, depending on the preparation protocol. 
At fixed a, the volume fraction at which a jammed sys- 



tem unjams and starts flowing (<fiu) is greater than the 
volume at which a flowing systems jams (<fij). A simple 
one-dimensional model, able to reproduce the observed 
features, clarifies that the hysteresis is due to the iner- 
tia of the system, and to the presence of an underlying 
energy landscape with many minima. 

At small by finite a, where the hysteresis is negligible, 
the jamming transition has a mixed first-order second- 
order character. Continuous quantities include P, G as 
well as the mean contact number Z (at variance from the 
(7 = case [lj]). Discontinuities are found, in close anal- 
ogy with thermal glass forming systems, in the asymp- 
totic value of the self-intermediate scattering function 
(the non-ergodicity parameter f^) and of others struc- 
tural relaxation functions. 

Numerical Model - We consider a system of spheres of 
diameter d and mass m in a box of size l x = l y = I6d 
and l z = 8d (we have investigated values of l z up to 
64d to check for absence of finite size effects). We use 
periodic boundary conditions along x and y, while along 
z the system is confined by rough plates, made by a set of 
glued particles. The bottom plate is fixed, while the top 
one (with mass ml x l y /d 2 ) is subject to a shear stress a. 
Particles interact via the linear spring dashpot model [|[ . 
When the center-center separation ry between particles 
i and j is smaller than d (S = d — > 0), particles i 
and j interact via a repulsive normal force 

F„ = k n Sn - 7„v„ (1) 

where n = r^/|r^-|, v n = ndfaj ■ n)/dt, and j n chosen 
fixing the restitution coefficient to e = 0.88 (the value of 
e does not qualitatively change our results). Particles are 
also subject to a viscous force — 77 so iv acting at all times, 
and different values of ?7 so i are explored. Lengths, masses, 
times, velocities and stresses are expressed in units of: d, 
m, t = y^m/kn, v — do/t, <tq = k n /d. Note that the 
shear stress is measured in units of the Young modulus 
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FIG. 1: (Color online) Dependence of the shear viscosity of 
the shear modulus on the volume fraction in the presence 
(panel a, high a) and in absece (panel b, small a) of hystere- 
sis. In the prence of hysteresis, the volume fraction at which 
the viscosity diverges is greater than that at which the shear 
modulus vanishes. Plain lines are power law fits. Panel (c): 
jamming phase diagram at T — 0. Full symbols mark the 
transition from the jammed to the flowing state (the shear 
modulus vanishes), while open symbols mark the transition 
from the flowing to the jammed state (the viscosity diverges). 
In the gray region, the system has an hysteretic behavior. 
The hysteresis decreases on increasing the viscosity 7 of the 
solvent the particles are immersed in. The inset illustrates 
that 4>j decreases as 7 increases, until it reaches <j>u- 



of the particles, implying that the hard spheres behavior 
is recovered in the a — > limit. 

The system is prepared via a Lubachevsky- 
Stillinger [ij [ll| like procedure. Particles are placed into 
the system with infinitesimal radii, rapidly increased 
to their final value. Then, the system is allowed to 
relax until the kinetic energy vanishes. As in Ref. [ill ], 
using periodic boundary conditions in all directions we 
find (/> rcp ~ 0.645 as the maximum volume fraction at 
which the system is able to relax in an unjammed (zero 
pressure) state fiol j . 

Hysteresis at jamming transition - We have found the 
response of the system not to be always determined by 
the volume fraction <j) and by the applied shear stress 
a. For some values of <j) and a, it also depends on the 
preparation protocol. We have therefore considered two 
limiting protocols, where a is fixed either coming from 
the jammed or from the flowing phase, respectively. In 



the first case, the shear stress is slowly increased starting 
from a = 0, until the desired value is obtained. In the 
other case, a is first fixed to a high value at which the 
system is seen to flow, and then decreased to its final 
value. Once the system reaches the steady state, we have 
measured the shear viscosity 77 = v s /al z , where v s is the 
mean velocity of the top plate, and the shear modulus 
G. The shear modulus is G = l z Sa/5L, where where 5L 
is the displacement of top plate of a system of volume 
fraction (f> jammed under the action of a shear stress a, 
occurring when a perturbing stress 8a is superimposed. 

Fig.s QJi,b show the dependence of the shear viscos- 
ity and of the shear modulus on the volume fraction, at 
different values of the shear stress. Each point is the av- 
erage over 20 independent runs, and errors are smaller 
that the symbol size. Measures taken coming from the 
flowing state, see the viscosity to increase with <fi, until 
the system jams at some value (f>j. At this point, the 
system is a solid with a finite value of G. Measures taken 
coming from the jammed phase, see G to decrease with 0, 
until the system becomes a liquid at some value <f>u < <j>j. 
At this point, the system starts flowing with a finite vis- 
cosity While at a small value of a (panel a) we found 
~ <f>j, and higher values of a (panel b) the transition 
is clearly hysteretic. Practically, we have identified <j>j 
and 4>u via power law fits of the data shown in Fig. [T^. 
The vicosity appears to alway diverge as (<fij(a) — 4>)~ s 
with s ~ 1.5, while G vanishes as (cf> — <?V (c)) 7 with 
7 ~ 0.50, in agreement with previous results at a = PJ. 

We have obtained data as those shown in Fig. [T^,b for 
different values of the applied shear stress a, and we have 
therefore identified the jamming <fij(o~) and unjamming 
4>u{&) transition lines, which are shown in the <fi-o~ plane 
of Fig.QJ. We have checked that these lines do not change 
if estimated at fixed <j) varying a. Both lines ends at 
(j) ~ rcp (J point) in the a — > limit. The width of the 
region where hysteresis continuously decreases with the 
applied shear stress. 

The hysteresis depends on the viscosity of the fluid the 
particles are immersed in ?7 so i. In particular, the unjam- 
ming transition 4>u{4>) does not depend on 77 so i, as one 
may have expected since 4>u is determined coming from 
the jammed phase, where ?7 so i cannot play a role. On the 
contrary, <fij depends on ?y so i: at a fixed value of the ap- 
plied shear stress cr, <pj decreases as ?7 so i increases, until it 
reaches 4>u and the hysteresis vanishes. The dependence 
of (fij and 4>u on ?y so i is shown in the inset of Fig. [T]:, for 
a = 0.17. 

These results suggest a physical interpretation of the 
hysteretic behavior. In the region of the control param- 
eters where hysteresis occurs, there are particle config- 
urations (energy minima) which are able to sustain the 
applied stress. If the system reaches one of these con- 
figurations coming from the jammed phase, it has zero 
kinetic energy, and stays jammed. Conversely, if the sys- 
tem reaches one of these minima coming from the flowing 
state, it has a kinetic energy large enough to escape from 
the minima, without jamming. 
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FIG. 2: (Color online) The same quantities shown in Fig. [T] 
are here shown as predicted by the model of Eq. [2] where / 
and V play the role of the applied stress a and of — <f>rc P , 
respectively. The model reproduces the hysteretic behavior, 
as well as the dependence of the width of the hysteretic region 
on both the applied stress and the viscosity of the solvent r) so \ . 



Model - The hysteretic behavior observed at the jam- 
ming transition, and its dependence on rj so i, are qualita- 
tively explained by a simple model for the motion of a 
particle in a energy landscape with many minima W(x), 
subject to a constant driving force / and to a viscous one 



mx = f — dW(x)/dx — r) so \X. 



(2) 



For simplicity sake, here we use W(x) = — V cos{ojx), 
and fix to = 1 and lu = 2tt. V measures the depth of the 
energy minima, and we therefore relate it to the volume 
fraction (V increases as cf> overcomes </> rcp ). With this 
choice for W(x), Eq. [2] has been extensively investigated 
as it models a large variety of physical processes (see 
Ref. [To| for a comprehensive review) . Eq. [2] admits two 
type of solutions, known as 'running' (the mean velocity 
is (x) > 0) and as 'locked' (x = 0), which correspond to 
our flowing and jammed state, respectively. In the run- 
ning state, it is possible to define the shear viscosity as 
i] = (x)/f. In the locked state, the equilibrium position 
x eq is fixed by / = dW(x)/dx\ x , and the shear mod- 
ulus is G(f,V) = d 2 W/dx 2 \ x<! . Eq. [2] reproduces the 
hysteretic behavior observed across the jamming transi- 
tion because in a given of the control parameters V and 
/ both the running and the locked solution are allowed, 



the observed one depending on the preparation protocol. 
For instance, we show in Fig.s [2Kb the dependence of 
the shear viscosity and of the shear modulus on V. The 
shear viscosity, measured when V increases, diverges at 
a jamming threshold Vj, while the shear modulus, mea- 
sured decreasing V, vanishes at an unjamming thresh- 
old Vu < Vj. This behavior closely resemble that of 
Fig. [Ha). 

From Eq. [2] it is also possible to determine the 'jam- 
ming diagram' in the V-f plane. The unjamming line 
is that where the shear modulus vanishes, G(f, vq) = 
(i.e. Vu = u/f), while the jamming line (numerically de- 
termined) is that where the shear viscosity diverges. The 
resulting diagram, which is shown in Fig. HJc), qualita- 
tively reproduces that of Fig. [IJc). In particular, the 
hysteretic region decreases as / decreases. Moreover, as 
shown in the inset, Eq.[2]also reproduces the dependence 
of the width of the hysteretic region on the shear viscos- 
ity. 

Structural signatures of the jamming transition - We 
now consider how dynamical and structural quantities of 
the system changes at the jamming transition. We fo- 
cus on small values of the applied shear stress a < 0.1 
where, as shown in Fig. [TJb), hysteretic effects are negli- 
gible. We show in Fig. [3Ja) the dependence of the mean 
contact number on the volume fraction: Z continuously 
increases with <fi, and has a cusp at the transition. Above 
the transition, the effect of the applied shear stress on Z 
is negligible, and Z grows as a power law with exponent 
~ 0.5 as at a = The pressure P, shown in Fig. [3^b) 
also increases with <j>, linearly above the jamming thresh- 
old. Note that, due to the continuous collisions of the 
flowing grains, below the jamming transition both the 
pressure and the mean contact number have finite val- 
ues, at variance with the a — case [l| . The a = case 
is approached continuously as a decreases. For instance, 
we have evaluated via a polynomial fitting the left deriva- 
tive of Z at 4>j, lim^^ dZ/d<j>, finding that it diverges 
as cr~ x , x ~ 0.2, as a decreases. 

While Z , G and P (not shown) are continuous at the 
transition, there are quantities which are discontinuous. 
Precisely, the asymptotic values of some two-time cor- 
relation functions are discontinuous, as in glass form- 
ing systems. Two of them are shown in Fig. [3Jd and 
Fig. Et- The first one is the self-intermediate scat- 
tering function F(k, t) — (3>(k, t)), where $(k, t) = 
1/N Y^j exp(— ik • (rj (t) — (0))) , we have investigated in 
the direction k = (0,27r/d, 0) perpendicular to both the 
confining plates and to the driving force. F(k, t) goes 
to zero in the flowing regime, while it stays to one when 
the system is jammed. The asymptotic value of F(k,t), 
known as non-ergodicity parameter , is therefore dis- 
continuous at the transition. The second quantity is the 
two-time mean contact number 



Z(t) 
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(3) 



where the sum runs over all particles, and Cjj(i) = 1(0) if 
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FIG. 3: (Color online) Structural changes at the jamming 
transition at small a. Panel (a) and (b): dependence of the 
mean contact number and of the normal pressre on the vol- 
ume fraction, for different values of a, as indicated. Panel 
(c): self-intermediate scattering function for different values 
of (f> (main panel), and volume fraction dependence of its 
asymptotic value /oo (inset). Panel (d): two-time mean con- 
tact number Z(t)/Z(0) (main panel) for different values of 
(j> (same as in panel c), and <f> dependence of its asymptotic 
value Z(oo)/Z(0) (inset). In panels a,c and d, a = 2 10 -3 . 



particles i and j touch (do not touch) at time t. Z(0) is 
the usual mean contact number. As shown in Fig. [3ji, 
Z(t)/Z(0) goes to zero when the system flows, while 
stays to one when the system is jammed. Its asymp- 
totic value is therefore discontinuous at the transition, 
exactly as /oo. The analogy with the glass transition 
is also validated by the study of the relaxation time r 



(F(k, t) — 1/e), which we have found to diverge with 
the shear viscosity r ~ ij ~ {<pj — (f>)~ s (s ~ 1.5) as 
the transition is approached, and by the study of the dy- 
namical susceptibility %4, defined as the fluctuation of 
the self-intermediate scattering function, having a maxi- 
mum x* at a time t* , with t* cx 77, and x* <x {4>J — < t } )~ u 1 
v ~ 0.8. Since our system relaxes at densities which are 



close to 1 



shear appear to be much more efficient than 



thermal motion in inducing the structural relaxation of 
the system. The scenario of Fig. [31 which we have found 
at small a and ?7 so i = 0, is also observed at high a pro- 
vided that 77 so i is large enough for the hysteresis to be 
negligible. 

Discussion - The hysteretic nature of the jamming 
transition in the 0-er plane implies the presence of mem- 
ory, not predicted in the frictionless case considered here. 
At small but finite a, where the hysteresis is negli- 
gible, the jamming transition has a mixed first-order 
second-order character which is close to that observed 
in glass forming systems. One time quantities, including 
the mean contact number, are continuous at the tran- 
sition, but the asymptotic value of two-time correlation 
functions, such as that of the self-intermediate scatter- 
ing function (the non-ergodicity parameter) or that of 
the two-time mean contact number, is discontinuous. 
The jamming transition at a = appears therefore the 
(7^0 limit of a glassy-like transition observed at all 
values of a. Important open questions ahead include 
the understanding of the role of frictional forces, which 
must be taken into account to describe real granular sys- 
tems [ill, 13 , jrl - EEI, [H| , as well as that of a finite tem- 
perature 171 Hsf, 



We thank A. Fierro and M. Nicodemi for interest- 
ing discussions, and acknowledge computer resources 
from CINECA, Unina-SCOPE Grid, INFN-Grid and 
CASPUR. 



[1] C.S. O'Hern et ai, Phys. Rev. Lett 88, 075707 (2002); 

C.S. O'Hern et al, Phys. Rev. E 68, 011306 (2003). 
[2] A. J. Liu, S. R. Nagel, Nature 396, 21 (1998). 
[3] J. Schwarz, A.J. Liu and L.Q. Chayes, Europhys. Lett. 

73, 560 (2006). 
[4] C.F. Moukarzel, Phys. Rev. Lett. 81, 1634 (1998). 
[5] J.N Roux, Phys. Rev. E, 61, 6802 (2000). 
[6] P. Olsson and S. Teitl, Phys. Rev. Lett. 99, 178001 

(2007). 

[7] C. Heussinger and J.L. Barrat, Phys. Rev. Lett. 102, 

218303 (2009). 
[8] L. E. Silbert et al, Phys. Rev. E 64, 051302 (2001).. 
[9] B.D. Lubachevsky and F.H. Stillinger, J. Stat. Phys. 60, 

561 (1990). 

[10] H. Risken, The Fokker-Plank Equation, Springer- Verlag, 
1984. 

[11] H.P. Zhang and H.A. Makse, Phys. Rev. E 72 011301 
(2007). 



[12] M. Pica Ciamarra, M. Nicodemi and A. Coniglio, Phys. 

Rev. E 75, 021303 (2007). 
[13] F. Lechenault, O. Dauchot, G. Biroli, J. -P. Bouchaud, 

EPL 83, 46003 (2008). 
[14] E. Somfai et al. Physical Review E 75, 020301 (2007). 
[15] D.S. Grebenkov et ai, Phys. Rev. Lett. 100, 07801, 

(2008). 

[16] M. Pica Ciamarra et al, in preparation. 

[17] Z. Zhang et ai, Nature 459, 230 (2009). 

[18] L. Berthier and T.A. Witten, arXiv:0903 

[19] A corrective term is introduced to take into account 
the finite size of our sample. In the presence of the 
rough confining plates, we define the volume fraction as 
cj>(N) = Nv/lJ y l z [1 + 5V/(lJ y l z )], where v is the vol- 
ume of a grain, and we have introduced a term SV/l x l y l z 
to take into account the effects of the walls protruding 
into the system. We have determined SV ~ 52 requiring 
the maximum volume fraction at which the system re- 



■5 



laxes in a zero pressure state in the presence of walls to 
be equal to that obtained without walls. 



